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Abstract 

We consider a Friedmann-Robertson-Walker spacetime filled with both viscous radiation and 
nonviscous dust. The former has a bulk viscosity which is proportional to an arbitrary power of 
the energy density, i.e. £ oc p”, and viscous pressure satisfying a nonlinear evolution equation. 
The analysis is carried out in the context of dynamical systems and the properties of solutions 
corresponding to the fixed points are discussed. For some ranges of the relevant parameter 
v we find that the trajectories in the phase space evolve from a FRW singularity towards an 
asymptotic de Sitter attractor, confirming and extending previous analysis in the literature. For 
other values of the parameter, instead, the behaviour differs from previous works. 


Introduction 

The ACDM model, which is essentially a homogeneous and isotropic Friedmann-Robertson-Walker 
(FRW) spacetime filled with both pressureless dust and a cosmological constant, fits very well the 
evolution of the observable Universe from structure formation to the present accelerated expan¬ 
sion. The possibility of deviating from such a simple yet effective model has been widely considered 
throughout the literature, e.g. by taking into account less symmetric backgrounds or by introduc¬ 
ing less ideal matter sources. For instance, in the context of generalizing the type of fluid that 
sources the cosmological evolution, one could think of introducing dissipative terms in the energy- 
momentum tensor, in order to take into account nonequilibrium effects in the fluid(s). A first 
approach for describing nonequilibrium thermodynamic effects in a relativistic context was given 
by Eckart |lj. Unfortunately such a theory presents noncausal features, admitting superluminal 
propagation of the dissipative signals: this pathology lead other authors to seek for a second-order, 
causal extension, which could be identified in the so-called Israel-Stewart’s theory (IS) [2] (different 
approaches to the problem have been put forward, e.g. ES], and have been shown to agree with 
IS not far from equilibrium). In the present paper we work along these lines, by including in the 
dynamics of an expanding metric a nonequilibrium contribution to the equilibrium pressure p of 
the fluid, in the form of a bulk viscous pressure II (some previous studies on the topic can be found 
e.g. in [6j [5, 7j [SJ [9, (TO) ITT, 12] [13], ll]). It is well known that, in considering a homogeneous 
and isotropic background, other dissipative processes with directional character - heat transfer and 
shear viscosity - do not contribute to the dynamics. The model governing the evolution of bulk 
viscous pressure that we choose to consider is a nonlinear extension of Israel-Stewart’s model (nIS 
from now on), presented for the first time in m and further analysed in ngnz]. Apart from the 
relaxational time r typical of IS causal theory, the nIS introduces a characteristic time r* for nonlin¬ 
ear effects such that, when r* —> 0, one recovers IS. The phenomenological nature of such a model 
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has been made clear in the original work and the general conditions underlying the hydrodynamic 
description (at the foundation of all the previously mentioned models) are also understood: for 
instance, one has to be aware that the requirement that the mean interaction time t c between fluid 
particles be less than the time scale 6~ 1 set by cosmological expansion (where 6 is the expansion 
scalar) could cease to hold in the case of accelerated expansion. 

The main task of the present work, carried out by means of a dynamical system approach, is 
to describe the space of solutions of an expanding FRW spacetime in presence of a double-fluid 
source: a pressureless dust component and a viscous radiation component whose bulk pressure is 
governed by the nIS model. We do not necessarily focus on the viability of the results in terms 
of observational consistency; rather we choose to explore in a more general way the impact of the 
specific nonlinear model on the past and future dynamics. 

We first present the background setting and the model in Sec{l] we then analyse the correspond¬ 
ing dynamical system in Secj2j where critical points of the system are identified and their stability 
assessed; in Secj3]we give interpretation of the results in terms of cosmological models and finally 
we give some concluding comment in Sec (4} 


1 The model 


Our analysis takes place in standard general relativity, hence Einstein’s field equations (EFE) 


G — n T) w 


are considered a valid starting point. From now on we will adopt the convention c = k = 1. 
The assumptions of homogeneity and isotropy justify the choice of a Friedmann-Robertson-Walker 
(FRW) background metric 

ds 2 = — dt 2 + a 2 (t ) ( dr 2 + r 2 dd 2 + r 2 sin 2 6 dcj ) 2 ) 

and specify in this way the form of the left hand side of EFE. Through the scale factor a(t) we 
can define the expansion scalar 9 = 3 d/a, where dot represents the derivative w.r.t. cosmic time 
t. The right hand side of EFE is built by choosing a noninteracting mixture of dust, with energy 
density pd, and viscous radiation, with energy density p v . Using the equation of state for radiation 
Pv = \pv, the latter fluid has an effective pressure p e ff = \p v + II, where the first term is the 
equilibrium part whereas II is the nonequilibrium contribution, i.e. the bulk viscous pressure. By 
specifying the metric and the matter content of the model we are then able to write down the 
relevant equations for the system. First of all, the Friedmann constraint equation coming from 
EFE, 

Pd + Pv = -9 2 , (1) 

allows us to disregard the evolution of the dust component and hence to focus only on the evolution 
for the energy density of the viscous fluid: by imposing energy-momentum conservation and making 
use of the constraint we have that 

Pv + . ( 2 ) 

Moreover, from EFE we obtain the Raychaudhuri equation 

e = -\e 2 - l -{ Pv + m) , ( 3 ) 
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where we made use again of eq. 0. 

The introduction of a viscous pressure begs for a specification of its nature and behaviour 
through some evolution equation. The model of evolution we will consider here has been proposed 
in m and analysed also in |16] for a single fluid and in P2| = 


rll = -Q 9 - n 


(l + n 



-i 



e + 


7 C 

T 

T C 

T 


(4) 


This is a nonlinear extension of IS equation. The nonlinear character of the nIS model is expressed 
by the term in round brackets and it is governed by a relaxational time r*: if such characteristic 
time r* —> 0, we recover IS theory. The nonlinear term itself is inversely proportional to the entropy 
production rate in the fluid and hence, for thermodynamic consistency, has to be positive. The 
elements that appear in eq. 0 are defined as follows: 

• bulk viscosity: ( = Co Pv, with Co > 0 and v > 0, 


• linear relaxational time: r = C/(7 c & Pv), 

• nonlinear relaxational time: t * = k 2 r, with real k, 

• temperature: T = Topi 7 1 ^ 7 , by assumption of barotropicity of the fluid. 


We recall that Cb is the dissipative part of the speed of sound, to be taken into account along 
with the adiabatic part c s in the total expression for the speed of sound V 2 = c 2 + c 2 . From the 
definition of adiabatic speed of sound and the condition V 2 < 1, it is easy to find the bound c 2 < 2/3. 


2 The Dynamical System 


In order to write the evolution equations in the form of an autonomous system, we define the 
following dimensionless variables: 


3 Pv 


6 2 



and Z 


Co 

3VQ1-2V- 


(5) 


We define also a new time variable r through dt/dr = 3/9, whose action on the variables will be 
denoted by a prime. Evolution towards increasing r is then related to the evolution of an expanding 
universe. By deriving the definitions in eqs. 0 w.r.t. r and making use of eqs. 0 © , we can write 
down the autonomous system: 


Q' = 

n' = 


z' = 


-(i-fi) 

— 4 cl 

,n 2 


q + 3n 
n 


i + 


3 z n 


v-\ 


11 + 


+ 3 : 


n 


n - 8 ) -(i-ft)n 


-Zl.-l 


n + 3 (n +1) 


3 k 2 




( 6 ) 

(7) 

( 8 ) 


Due to the constraint eq.Q we have D E (0,1] and n E (—0 3c^/4A; 2 , + 00 ) from the requirement 
of positive entropy production. It is immediately clear from eq.Q that the system is ill-defined in 
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the subspaces Z = 0 and Cl = 0 so, apart from locating the critical points, we will need to study 
the behaviour of the system in a neighborhood of such planes in the phase space: this will be done 
through a numerical analysis. Note also that the singular behaviour in Z = 0 will prevent this 
variable to change sign during the evolution: in the following we will focus on the part where Z 
> 0. 

In the case when v = 1/2, eq.Q gives Z = constant, as can be seen also from the definition in 
eq.([5]): the dimensionality of the system is thus reduced from 3 to 2 and the analysis is formally 
equivalent to the one carried out in m- With regard to this case we will just report some general 
considerations in Sec. 12.31 


2.1 Finite analysis (v ^ 1/2) 

We start by locating and analyzing the critical points with finite values of the coordinates in the 
phase space. 


As long as v ^ 1/2, we note that a generic finite critical point {f2 c , IT C , Z c } with Z c ^ 0, if it 
exists, describes a de Sitter model because, from the definition of Z, 


6 = 



(9) 


and constant. As a consequence the scale factor of the model in that point will be oc Exp ($ot/3). 
The energy density and viscous pressure in such a model would attain, respectively, the constant 
values 


1 

3 


ci c 0 q 


and 


n c &l 


( 10 ) 


In the invariant set II = 1 (i.e. where Cl' = 0) we find a finite critical point of this kind: 


Pn = 



4 cl 


9 (<£-1/2) (1 ~k 2 /c 2 ) 


( 11 ) 


Requiring positivity of Z imposes either {c^ > 1/2 A k 2 < c^} or {c 2 < 1/2 A k 2 > c 2 }, but only 
in the former case Po lies in the physical region of the phase space given by the constraints on the 
variables. Thus, considering the bound on Cf, found previously, the existence of this critical point 
is subject to the condition 

1/2 < c 2 b < 2/3 . (12) 

Also in m possibility of exponential inflation was found in this range of velocities. The location of 
Po does not depend on u, but the stability properties depend crucially on the exponent. By study¬ 
ing the linearized system around the point we conclude that Pq is a stable attractor for v < 1/2 


and an unstable saddle for v > 1/2. In the following we will assume validity of the bound eq. (12), 
because the dynamics results to be richer. We will make some comments on the case 0 < c 2 < 1/2 
in Sec. [H 


A second finite critical point can be found by noting that one can solve the equilibrium condi¬ 
tion for the system by approaching II = Z = 0 in the invariant set Cl = 1 along the line II = — 3 Z. 
In such a way P\ = {1,0,0} results to be a critical point of the system. The point Pi is a saddle, 
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Point: 

Po 

Pi 

P 2 

P+ 

P~ 

0 < v < 1/2 

sink 

saddle 

saddle 

(n/a) 

(n/a) 

v = 1/2 

(n/a) 

(n/a) 

saddle 

source 

sink 

v > 1/2 

saddle 

saddle 

sink 

(n/a) 

(n/a) 


Table 1: Stability of the finite critical points of the system. 


so it will act only as a transient stage for some trajectories in the phase space. 

It is also possible to discover another critical point in the finite phase space if we approach the 
line II = Z = 0 along a generic curve II oc Z n (with n > 1) and then we let II —> 0. From numerical 
analysis, the point P -2 = { 0 , 0 , 0 } found acts either as a saddle (if 0 < u < 1 / 2 ) or as a sink (if 
v > 1 / 2 ). 

The stability properties of the finite critical points are summarized in Table 0 In Fig{l] some 
representative trajectories in the phase space are plotted, highlighting in this way the impact of the 
choice of v on the overall behaviour. The shaded part of the phase space is the negative entropy 
production region and trajectories emerging from it are denoted by dotted lines. The direction of 
the trajectories is specified by arrowheads and the locations of the finite critical points are identified 
by black dots. From the plots (and from the absence of sources for v 7 ^ 1/2 in Table 0) it is clear 
that the asymptotic points play an important role, as we will see in the next section. 


2.2 Asymptotic analysis (u 7 ^ 1/2) 

The study of the critical points at infinity is carried out by compactifying the phase space. This is 
done by redefining the unbounded variables through the following transformation: 

f 7 * 

II =- cos0 , Z =- siiuj) , (13) 

1 — r 1 — r 


where r E [0,1) and 4> £ [0, 7 r]. In this way, the directions of divergence of (either one or both) the 
variables are mapped on the boundary r = 1. Note that, because of the lower bounds on II and 
Z, at infinity the variable 4> is restricted to [0, 7t/2] . By using eq.(13) we obtain a system in the 
implicit form 


n' = f(n,r,cj>) 

r’ = 9 ( 12 , r, <j>) 

</> = h ( 12 , r, <j>) . 

The location of the critical points at infinity can be obtained by evaluating first the limit r —> 1~ 
and then by finding the values {O c , <f > c } that satisfy Q' = (/>' = 0. The stability of the critical points 
is then checked through the linearization matrix of the 2 -dimensional phase space { 12 , 0 } at infinity 
together with the sign of r\ which reflects the stability in the radial direction. We made use of 
numerical studies in order to assess the stability of some of the asymptotic critical elements. 
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Figure 1: Representative trajectories in the phase space for q, = 0.8 and k = 0.4. Arrowheads 
specify the direction of the trajectories. The location of Pq is specified by a black dot: for u = 0 
(left) it is a sink, for v = 3/2 (right) it is a saddle. The region of negative entropy production is 
shaded and trajectories emerging from it are dotted. 
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Point: 

poo 

Lf 

700 

l 2 

777 , 2 ° 

0 < v < 1/8 

source 

(n/a) 

saddle 

saddle 

1 / = 1/8 

source 

saddle 

saddle 

saddle 

1/8 < v < 1/2 

source 

(n/a) 

saddle 

saddle 

v = 1/2 

(n/a) 

(n/a) 

(n/a) 

(n/a) 

1 / > 1/2 

saddle 

(n/a) 

source/sink 

saddle 


Table 2 : Stability of the asymptotic critical elements of the system. 


For r —> 1 , the system at infinity is given by 

14 7 = — 3cos <j>- -- (14) 

1 — r 

/ 3 cos <f> 

8V~ ' 

■ ^5 cos 2 <fi + 4 — - — 2 fu + cos(214^ (15) 

= sni^-r) wsin0 ( 5 ~ 8 + I) n ) ' (16) 

First of all, in the particular case when v = 1/8 we note that the condition Q' = cf) 1 = 0 
identifies the value 14 = 1, V</>: thus for this choice of the exponent we find an asymptotic critical 
line LJ° = { 1 ,(/)} with saddle behaviour. 

For v 7 ^ 1/8 there is no critical line L but a critical point in P£° = {1,0}: from numerical 
studies this point turns out to be a source for v < 1/2 and a saddle for v > 1 / 2 . 


Finally, the common factor cost/ in eqs.(14) and (16) tells us that the system has a critical line 

= {14,7r/2}, which is a saddle for 14 E (0,1). In this line we can identify two particular critical 
points: 1 3 ° = {1,7r/2} and rrt^° = { 0 , 7 t/ 2 }. 

The point is equivalent to Z — >• 00 in the invariant set 14 = 1. By solving II' = 0 for Z 
with 14 = 1 in the original system, it is possible to see that this asymptotic critical point actually 
corresponds to two vertical asymptotes II = ±^\/2 |c&| in the invariant set, with different stability 
features: along the positive asymptote, is a saddle for v < 1/2 and a source for v > 1 / 2 ; along 
the negative asymptote, it is a saddle for v < 1/2 and a sink for v > 1/2. The degeneracy in 
this asymptotic point is made apparent in the bottom panel of Fig j3j representing the compactihed 
invariant set: the point’s double nature of source and sink is shown for a representative trajectory 
in a case where v = 3/2. 

The point m 2 ° corresponds to Z — > 00 in the plane 14 = 0. This point is a saddle Vza 


2.3 The case v = 1/2 

If the exponent in the definition of the bulk viscosity takes the value 1/2, we have a drastic 
modification of the system: in fact, in this case the variable Z becomes a constant and the phase 
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space is subject to a reduction in dimensionality. The dynamics now takes place in 2-dimensional 
invariant subsets (Z 7 = 0) of the original 3-dinrensional problem, one subset for each value of Z. 
Defining the constant Zq = Co/3 1 / 2 , we have the following system of equations: 


n' = 

n' = 


o + 3n 
fm 1 / 2 


-(l-D) 

—4 c 2 n 

+ 3^ (n - ^ ) - (i-n)n 


i + 


3 Zq 


n + 


3 k 2 

4 4 



(17) 


(18) 


The situation is analogous to the case treated in HZ], the only difference being the initial choice 
for the bulk viscosity C given in Sec. [l] Nevertheless the system presents three critical points, 
corresponding to three roots of a cubic equation: 

• P + is a source 


• P is a sink 

• P 2 is a saddle 

Both the source and the sink represent polynomially expanding solutions with viscous fluid pre¬ 
dominance. A future de Sitter attractor is possible but a fine tuning of the parameters is needed. 
The details in |17| are qualitatively equivalent to this case and we will not discuss them further. 


3 Cosmological evolutions 

In translating the properties of the phase space elements in terms of cosmological models we choose 
to disregard the trajectories emerging from the boundary of the unphysical region (shaded in gray 
in all the Figs.). The main feature of their past stage is obviously the fact that the entropy pro¬ 
duction rate diverges and the system is not well defined. Focusing on the rest of the phase space, 
we analyse the cosmological models corresponding to different ranges of v. 

When 0 < v < 1/2, the main elements are a source in P^° and a sink in Pq. Thus, for a —> 0 
we have a radiation-dominated FRW model, while for a —> 00 the system evolves towards a stable 
de Sitter model. Possible intermediate stages are a nonviscous dust dominated era in P 2 and a vis¬ 
cous radiation dominated era in l ™, both characterized by polynomial expansion. Such transient 
stages are identifiable through the behaviour of the solid trajectories in the left panel of Fig|l] The 
compactified invariant set ST = 1 is represented in the top panel of Fig. [2] When crossing the value 
v = 1/8, there is a change in the relative attraction between P^° and 1%?: if v < 1/8 the direction 
of the trajectories at infinity is P^° —> , while for v > 1/8 is the opposite. In the latter case, as a 

consequence of this change, we also note that the point Pf° (apart from being a source) can act as 
a sink for trajectories lying in the invariant set, as shown in the bottom panel of Fig. [2] This means 
that there exist solutions that are FRW both for a —> 0 and a —> 00, but this is obviously an unsta¬ 
ble particular situation if the fluid is not purely viscous, because of the instability in the fl direction. 

If v > 1/2 the stability of the critical points changes drastically: the past attractor for the tra- 

1 

jectories is now 1 2 ° along the positive asymptote, with a scale factor proportional to (t — to) 2 ( 1 +'/2 l%l) 
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Figure 2: Compactified phase space in the invariant set 0 = 1 with c b = 0.8, k = 0.4 and v = 0 
(top) or i/ = 9/40 (bottom). The critical points are denoted by black dots. Negative entropy 
production subspace is shaded in gray, while the region of stability in the O direction is shaded in 
blue. 
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Figure 3: Compactified phase space in the invariant set 0 = 1 with q, = 0.8, k = 0.4 and v = 3/2 
(top). The critical points are denoted by black dots. Negative entropy production subspace is 
shaded in gray, while the region of stability in the 11 direction is shaded in blue. Detail of the 
region around l $° for the same values of the parameters (bottom): the representative trajectory 
shows the double nature (source/sink) of the critical point. 
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which denotes a slower approach to a —> 0 with respect to the pure radiation case. The main sink 

for a —> oo is P 2 , which is a polynomially expanding model dominated by non-viscous dust with 

a oc (t — f 0 ) 2/3 . As we have seen, the point 1 2 ° itself is also a sink for some trajectories along the 

negative asymptote. The sink in is characterized by a —> 00 , p v —> 00 and |II| —> 00 , which 

_ 1 

could represent a future Type I singularity (or Big Rip [18]) with scale factor oc ( t. s — t ) . 

We can calculate the deceleration q and effective EoS 7 e ff parameters in this point, finding 


q = 1 - 2 v / 2 |c fo | 

7 eff = o (f - V2\Cb 


(19) 

( 20 ) 


Given the bound on Cb in eq.(12) we have that q < — 1 and 7 e // < 0, showing an effective 
phantom behaviour. As regards the de Sitter model Pq, being a saddle in this range of v, it acts 
as a possible transient stage. The situation in the compactified invariant set 17 = 1 is represented 
in the top panel of Fig. [ 3 ] for a value v = 3/2. 


Some considerations on the bound eq.(12) are in order. When —> 1/2, the de Sitter critical 
point is pushed towards infinity along the negative asymptote we discussed earlier, coinciding with 
^ 2 ° in the limit. In the range 0 < < 1/2, the de Sitter critical point Pq is effectively no longer part 

of our phase space and the point l “ along the negative asymptote inherits its stability properties: it 
is a sink for 0 < v < 1/2 and a saddle for v > 1/2. However, the different bound on q, changes also 


the cosmological nature of the point, as can be seen from eqs.(19)(20): we don’t find a phantom 
model because now — 1 < q < 1 and 4/3 < 7 e // < 0. 


4 Conclusions 


In the present paper we analysed the phase space of a FRW spacetime filled with both radiation 
and dust. The effective pressure of the radiative component splits in an equilibrium and a nonequi¬ 
librium part: the latter is a viscous pressure term which satisfies a nonlinear evolution equation, 
given by eq. 0- This study is a generalization of a previous work m to the case of a general 
power-law dependence of the bulk viscosity on the energy density, i.e. ( = Co Pv with v > 0. We fo¬ 
cused mainly on the case where the bound on velocities eq.( 12) is satisfied because in such a range a 
de Sitter critical point, interesting from the cosmological point of view, is present in the phase space. 


The system displays interesting features regarding the stability of the critical points (summa¬ 
rized in Tables [T] and [ 2 ]) and hence the evolution of the trajectories in the phase space. As already 
known in the literature, the stability properties of the system change when crossing the value 
v = 1/2. Our findings for 0 < v < 1/2 are in line with previous studies, such as j5j 6f (1(7. IT] : 
the evolution for an expanding universe starts with a FRW phase for a —> 0 and is asymptotically 
de Sitter for a —> 00 , a behaviour that survives in this model since the Eckart approach; in our 
double-fluid system also an intermediate matter dominated stage is possible. The overall evolution 
can be schematically represented in the following diagram: 
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The arrows point in direction of increasing scale factor and the sign in round brackets specifies the 
asymptote (see Sec. 2.2). If v > 1/2, differently from other works, we find a viscous radiation- 
dominated FRW behaviour for a —> 0; only in some intermediate stage can we find a de Sitter 
phase; finally, there are two possible futures for a —> oo, either a dust-dominated FRW model or a 
viscosity-driven phantom model: 
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